function hfile = sdm_AddFilters(hfile, opts)
% SDM::AddFilters  - add filtering regressors to SDM
%
% FORMAT:       [sdm =] sdm.AddFilters(opts)
%
% Input fields:
%
%       opts        settings
%        .constant  also add constant regressor (default: false)
%        .ftype     either of 'dct', {'fourier'}, 'linear'
%        .number    number of filters/frequencies (cycles, 3)
%        .timepts   number of time points (valid only for empty SDM)
%
% Output fields:
%
%       sdm         altered SDM with added regressors

% Version:  v0.8a
% Build:    9102000
% Date:     Oct-20 2009, 12:46 AM CEST
% Author:   Jochen Weber, SCAN Unit, Columbia University, NYC, NY, USA
% URL/Info: http://wiki.brainvoyager.com/BVQXtools

% argument check
if nargin < 1 || ...
    numel(hfile) ~= 1 || ...
   ~isBVQXfile(hfile, 'sdm')
    error( ...
        'BVQXfile:BadArgument', ...
        'Invalid call to %s.', ...
        mfilename ...
    );
end

% get design matrix right
bc = bvqxfile_getcont(hfile.L);
sdm = bc.SDMMatrix;

% options check
if nargin < 2 || ...
   ~isstruct(opts) || ...
    numel(opts) ~= 1
    opts = struct;
end
if ~isfield(opts, 'constant') || ...
   ~islogical(opts.constant) || ...
    numel(opts.constant) ~= 1
    opts.constant = false;
end
if ~isfield(opts, 'ftype') || ...
   ~ischar(opts.ftype) || ...
   ~any(strcmpi(opts.ftype(:)', {'dct', 'fourier', 'linear'}))
    opts.ftype = 'fourier';
else
    opts.ftype = lower(opts.ftype(:)');
end
if ~isempty(sdm) || ...
   ~isfield(opts, 'timepts') || ...
   ~isa(opts.timepts, 'double') || ...
    numel(opts.timepts) ~= 1 || ...
    isinf(opts.timepts) || ...
    isnan(opts.timepts) || ...
    opts.timepts < 4
    if isempty(sdm)
        error( ...
            'BVQXfile:BadArgument', ...
            'Empty SDM: Number of timepoints required.' ...
        );
    end
    opts.timepts = size(sdm, 1);
else
    opts.timepts = fix(opts.timepts);
    sdm = zeros(opts.timepts, 0);
    bc.PredictorNames(:) = [];
    bc.PredictorColors = zeros(0, 3);
end
bc.NrOfDataPoints = opts.timepts;
if ~isfield(opts, 'number') || ...
   ~isa(opts.number, 'double') || ...
    numel(opts.number) ~= length(opts.number) || ...
    any(isinf(opts.number) | isnan(opts.number) | opts.number < 1 | ...
     opts.number > floor(0.5 * (opts.timepts - size(sdm, 2))))
    opts.number = 3;
end
if numel(opts.number) == 1
    opts.number = 1:opts.number;
else
    opts.number = unique(opts.number);
end
opts.number = opts.number(:)';
nv = opts.timepts;

% temporal filter type
switch (opts.ftype)
    
    % DCT-based filtering
    case {'dct'}
        
        % build X
        X = zeros(nv, numel(opts.number));
        Xn = cell(1, numel(opts.number));
        Xc = zeros(numel(opts.number), 3);
        n = 0:(nv - 1);
        for dc = 1:numel(opts.number)
            X(:, dc) = cos(pi * (2 * n + 1) * opts.number(dc) / (2 * nv));
            Xn{dc} = sprintf('DCF-f%d', opts.number(dc));
            Xc(dc, :) = max(0, [240, 224, 224] - 3 * opts.number(dc));
        end
        
    % sin/cos set filtering
    case {'fourier'}
        
        % build X
        X = zeros(nv, 2 * numel(opts.number));
        Xn = cell(1, 2 * numel(opts.number));
        Xc = zeros(2 * numel(opts.number), 3);
        n = 0:(nv - 1);
        for pc = 1:numel(opts.number)
            X(:, 2 * pc - 1) = sin(pi * (2 * n + 1) * opts.number(pc) / nv);
            X(:, 2 * pc) = cos(pi * (2 * n + 1) * opts.number(pc) / nv);
            Xn{2 * pc - 1} = sprintf('sin-f%d', opts.number(pc));
            Xn{2 * pc} = sprintf('cos-f%d', opts.number(pc));
            Xc(2 * pc - 1, :) = max(0, [224, 240, 224] - 2 * opts.number(pc));
            Xc(2 * pc, :) = max(0, [224, 224, 240] - 2 * opts.number(pc));
        end
        
    % only detrend
    case {'linear'}
        
        % build X
        X = ztrans(1:nv);
        Xn = {'linear_trend'};
        Xc = [224, 224, 224];
end

% add constant too?
if opts.constant
    X(:, end + 1) = 1;
    Xn{end + 1} = 'Constant';
    Xc(end + 1, :) = [255, 255, 255];
end

% add to content
bc.FirstConfoundPredictor = min(bc.FirstConfoundPredictor, size(sdm, 2) + 1);
bc.SDMMatrix = [sdm, X];
bc.NrOfPredictors = size(bc.SDMMatrix, 2);
bc.PredictorNames = [bc.PredictorNames(:)', Xn(:)'];
bc.PredictorColors = [bc.PredictorColors; Xc];
bc.RTCMatrix = bc.SDMMatrix;
cp = find(any(bc.SDMMatrix ~= 0, 1) & all(diff(bc.SDMMatrix) == 0, 1));
if ~isempty(cp)
    bc.RTCMatrix(:, cp) = [];
    bc.IncludesConstant = 1;
    if numel(cp) > 1
        cp(end) = [];
        bc.SDMMatrix(:, cp) = [];
        bc.PredictorNames(cp) = [];
        bc.PredictorColors(cp, :) = [];
    end
else
    bc.IncludesConstant = 0;
end

% set back
bvqxfile_setcont(hfile.L, bc);
